/* *****************************************************************************
*
* Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
* All Rights Reserved.
*
* This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
* the contents of this file may not be disclosed to third parties, copied or
* duplicated in any form, in whole or in part, without the prior written
* permission of Silicon Graphics, Inc.
*
* RESTRICTED RIGHTS LEGEND:
* Use, duplication or disclosure by the Government is subject to restrictions
* as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
* and Computer Software clause at DFARS 252.227-7013, and/or in similar or
* successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
* rights reserved under the Copyright Laws of the United States.
*
***************************************************************************** */
      subroutine matnoise(a,lda,nlin,ncol)
      real a(lda,1)
c
      init = 1325
c$doacross
      do 35 i = 1,nlin
          a(i,1) = 0.0
   35 continue
      do 30 j = 1,ncol
         do 20 i = 1,nlin
            init = mod(3125*init,65536)
            a(i,j) = a(i,j) + (init - 32768.0)/16384.0
   20    continue
   30    continue
      return
      end
      
      subroutine matpower( nlin,ncol,mat,power,lda)
      real mat(lda,*), power(lda,*)
	real re, im

	    il = 1
	    power(il,1) = mat(il,1) * mat(il,1)
	    do j = 2, (ncol-1)/2
		re = mat(il,2*j-1)
		im = mat(il,2*j)
		power(il,j) = re * re + im * im
		power(il,ncol-j+1) = power(il,j)
	    end do
	    if( mod(ncol,2) .eq. 2) 
     $		power(il,ncol) = mat(il,ncol) * mat(il,ncol)
c$doacross local(i,j,re,im)
	do j = 2, ncol
	    do i = 2, (nlin-1)/2
		re = mat(2*i-2,j)
		im = mat(2*i-1,j)
		power(i,j) = re * re + im * im
		power(nlin-i+1,j) = power(i,j)
	    end do
	end do
	    il = nlin
	    power(il,1) = mat(il,1) * mat(il,1)
	    do j = 2, (ncol-1)/2
		re = mat(il,2*j-1)
		im = mat(il,2*j)
		power(il,j) = re * re + im * im
		power(il,ncol-j+1) = power(il,j)
	    end do
	    if( mod(ncol,2) .eq. 2) 
     $		power(il,ncol) = mat(il,ncol) * mat(il,ncol)
	if ( mod(nlin ,2) .eq. 0) then
	end if
      return
      end

      subroutine matlog( nlin, ncol, power, lda)
      real power(lda,*)
	real t
c$doacross local(i,j,t)
	do j = 1, ncol
	    do i = 1, nlin
		t = power(i,j) + 1
		power(i,j) = LOG(t)
	    end do
	end do
      return
      end
	    
      subroutine matcolor( nlin,ncol,mat,color,lda)
      real mat(lda,*)
      integer color(ncol,*)
	real t, a, coeff,xx,yy
	integer red, green, blue

	xx = 0.0
	yy = 0.0
	do j = 1, ncol
	    do i = 1, nlin
		a = mat(i,j)
		if ( a .lt. xx) then
		    xx = a
		else if(  a .gt. yy) then
		    yy   = a
		end if
	    end do
	end do
	coeff = MAX(ABS(xx),ABS(yy))

c	print *, ' MIN =', xx, '   MAX =', yy

	if ( coeff .gt. 0.) coeff = .99 * 255. / coeff

c$doacross local(t, a, red,blue,green,i,j)
	do j = 1, ncol
	    do i = 1, nlin
		t= mat(i,j)
		a = ABS(t)
		red = coeff* (a - t) 
		if( red .gt. 255) red = 255
		green = coeff* a
		blue = coeff* (a + t)
		if( blue .gt. 255) blue = 255
		color(j,i) = (red) + (256*green) + (256*256*blue)
	    end do
	end do
      return
      end
	    
